# -*- coding: utf-8 -*-
"""
Created on 03 08 下午4:55 2022 
@Author : HHQUAN
@Email : 1075960398@qq.com

"""
# 自适应噪声抵消示例

import numpy as np
import matplotlib.pyplot as plt
from scipy import signal
import soundfile as sf

# 参考端语音信号
wav, fs = sf.read('/home/cl/PycharmProjects/Mytest01/方位估计/man_woman16k.wav')
plt.figure()
plt.specgram(wav+0.01*np.random.randn(wav.size), NFFT=960, Fs=16000, noverlap=480)
plt.colorbar()
# plt.show()
# 读取噪声信号
noise, samplerate = sf.read('Machinery114_16k.wav')

# 构建相关噪声信号
h = signal.butter(2, Wn=[2000, 4000], btype='bandpass', output='ba', fs=16000)
b, a = h[0], h[1]
noise2 = signal.lfilter(b, a, noise)
plt.figure()
plt.subplot(211)
plt.specgram(noise2, NFFT=960, Fs=16000, noverlap=480)
plt.subplot(212)
plt.specgram(noise, NFFT=960, Fs=16000, noverlap=480)
# plt.show()

# 参考信号
N = np.min(np.array([wav.size, noise.size]))
d = wav[0:N] + noise[0:N]
plt.figure(figsize=(12, 4))
plt.subplot(211)
plt.plot(np.arange(d.size)/fs, d)
plt.xlim(0, d.size/fs)
plt.subplot(212)
plt.specgram(d, NFFT=960, Fs=16000, noverlap=480)
# plt.show()

# 进行自适应噪声抵消
M = 16
x = np.zeros((M, ))
w = np.zeros((M, ))
mu = 0.4
d_hat = np.zeros((N, ))
e = np.zeros((N, ))
for n in range(N):
    # 细节1
    x[1:M] = x[0:M-1]
    x[0] = noise[n]
    y = signal.lfilter(b,a, x)
    d_hat[n] = np.dot(w, y)
    e[n] = d[n] - d_hat[n]
    w += mu*e[n]*y

plt.figure(figsize=(12, 4))
plt.subplot(211)
plt.plot(np.arange(e.size)/fs, e)
plt.xlim(0, e.size/fs)
plt.subplot(212)
plt.specgram(e, NFFT=960, Fs=16000, noverlap=480)

plt.show()

